//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-

/////////////////////////////////////////////////////////////////////////////
// declaration of Lapack routines 
////////////////////////////////////////////////////////////////////////////
#include "EigSolver/Lapackf.h"
#include <iostream>
using std::ios;

#if defined(ADD_)
#define dsyev dsyev_
#define ssyev ssyev_
#define dsyevx dsyevx_
#define ssyevx ssyevx_
#define dsyevd dsyevd_
#define ssyevd ssyevd_
#define dsygv dsygv_
#define zhegv zhegv_
#define ssygv ssygv_
#define zheev zheev_
#define cheev cheev_
#define dsytrd dsytrd_
#define dstebz dstebz_
#define dstein dstein_
#define dstegr dstegr_
#define dstevr dstevr_
#define dpotri dpotri_                                  
#endif

extern "C" {
  void dsyevx(const char* job, const char* range, const char* uplo,
               int* n, double* a, int *lda, double* vmin, double* vmax,
               int* imin, int* imax, double* tol,
               int* m, double* eigval, double* eigvec,
               int* ldz, double* work, int* lwork, int* iwork,
               int* ifail, int *info);
  void ssyevx(const char* job, const char* range, const char* uplo,
               int* n, float* a, int *lda, float* vmin, float* vmax,
               int* imin, int* imax, float* tol,
               int* m, float* eigval, float* eigvec,
               int* ldz, float* work, int* lwork, int* iwork,
               int* ifail, int *info);

  void dsyev(const char*, const char*, const int*,
              double*, const int*, double*, double*, int*, int*);
  void ssyev(const char*, const char*, const int*,
              float*, const int*, float*, float*, int*, int*);

  void dsygv(const int*, const char*, const char*, 
             const int* , double* , const int* , double* , const int* , double* , 
             double* , int* , int*);
  void zhegv(const int*, const char*, const char*, 
             const int* , complex<double>* , const int* , complex<double>* , const int* , double* , 
             complex<double>* , int* , double*, int*);
  void ssygv(const int*, const char*, const char*, 
             const int* , float* , const int* , float* , const int* , float* , 
             float* , int* , int*);

  void dsyevd(const char* job, const char* uplo, int* n,
              double* a, int* lda, double* w, double* work, int* lwork, 
              int* iwork, int* liwork, int* info);
  void ssyevd(const char* job, const char* uplo, int* n,
              float* a, int* lda, float* w, float* work, int* lwork, 
              int* iwork, int* liwork, int* info);

  void zheev(const char* job, const char* uplo, const int* n,
             complex<double>* a, const int* lda, double* w, 
             complex<double>* work, int* lwork, double* rwork, int* info);

  void cheev(const char* job, const char* uplo, int* n,
             complex<float>* a, int* lda, float* w, 
             complex<float>* work, int* lwork, float* rwork, int* info);

  void dsytrd(const char* , int* , double*, int*, double*, double*, 
              double*, double*, int*, int*);

  void dstebz(const char* range, const char* order, int* n, double* vl, 
              double* vh, int* il, int* iu, const double* tol, double* d,
              double* e, int* m, int* nsplit, double* w, int* ib,
              int *is, double* work, int* iwork, int* info);

  void dstein(int* n, double* d, double* e, int* m, double* w,
             int* ib, int* isp, double* z, int* ldz, double* work, 
             int* iwork, int* ifail, int* info);
  void dstegr(const char* , const char*, int *, double*, double* e,
              double* vl, double* vu, int *il, int* iu, const double* tol,
              int *m, double* w, double* z, int *ldz, int* isuppz, 
              double*, int *, int *, int *liwork, int *info);
  void dstevr(const char*, const char*, int *n, double* d, double* e,
              double* VL, double* VU, int *IL, int *IU, const double* ABSTOL,
              int * M, double* W, double* Z, int *LDZ, int *ISUPPZ, 
              double* WORK, int *LWORK, int *IWORK, int *LIWORK, int *INFO);
  void dpotri(const char*, int*, double*, int*, int*);
}

void 
OhmmsNumerics::diag(const char jobz, const char uplo, const int n, double* a, 
          const int m, double* eval, double* workd, int lwork,  int& info) {
  dsyev(&jobz, &uplo, &n, a, &m, eval, workd, &lwork, &info);
}

void 
OhmmsNumerics::diag(const char jobz, const char uplo, const int n, complex<double>* a,
          const int m, double* eval, complex<double>* workd,
          int lwork,  double* rwork, int& info){
  zheev(&jobz, &uplo, &n, a, &m, eval, workd, &lwork, rwork, &info);
}


void 
OhmmsNumerics::diag(const char jobz, const char uplo, const int n, float* a, 
          const int m, float* eval, float* workd, int lwork,  int& info) {
  ssyev(&jobz, &uplo, &n, a, &m, eval, workd, &lwork, &info);
}

void 
OhmmsNumerics::diag(const int itype, const char jobz, const char uplo, 
          const int n, double* a, 
          const int m, double* b, const int lm, double* eval, 
          double* workd, int lwork,  int& info) {

  dsygv(&itype, &jobz, &uplo, &n, a, &m, b, &lm, eval, workd, &lwork, &info);

  //for(int i=0; i<m; i++)
  //cout << eval[i] << endl;
}


void 
OhmmsNumerics::diag(const int itype, const char jobz, const char uplo,
                    const int n, complex<double>* a,
                    const int m, complex<double>* b, const int lm, 
                    double* eval, complex<double>* workd, int lwork,  double* rwork, 
                    int& info) {
  zhegv(&itype, &jobz, &uplo, &n, a, &m, b, &lm, eval, workd, &lwork, rwork,&info);

}

void 
OhmmsNumerics::diag(const int itype, const char jobz, const char uplo, 
          const int n, float* a, 
          const int m, float* b, const int lm, float* eval, 
          float* workd, int lwork,  int& info) {
  ssygv(&itype, &jobz, &uplo, &n, a, &m, b, &lm, eval, workd, &lwork, &info);
}
/***************************************************************************
 * $RCSfile: Lapackf.cpp,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: Lapackf.cpp,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

